!> author: 左志华
!> date: 2022-08-17
!>
!> Navier-Stokes Solver
!> Navier-Stokes 求解器
module sph_navier_stokes_solver

    use sph_kinds, only: rk
    use sph_region, only: region_t, update_virtual_particles
    use sph_linked_list_nnps, only: build => build_grid_from_loc, query, estimated_capacity
    use sph_list_type, only: list_t
    use twin_sph_list_type, only: twin_list_t
    use sph_density_summation, only: density
    use sph_smoothed_kernel_function, only: skf
    use sph_internal_force, only: internal_force
    use sph_external_force, only: external_force
    use sph_artifical_force, only: artificial_viscosity, art_visc
    use sph_env, only: hsml, nums
    use sph_command_line, only: cli

    ! 项目常量，广泛用于各个例程访问
    integer :: nsave                !! 保存步数
    real(rk) :: dt                  !! 时间步长

    abstract interface
        !> Navier-Stokes equation <br>
        !> 用于处理由当前粒子域状态求解力 (加速度)，以驱动粒子域运动和状态更新 (NS方程)
        subroutine NS_fcn(region, istep, t)
            import :: region_t, rk
            type(region_t), intent(inout) :: region     !! region of particles <br>
                                                        !! 粒子宏观变量域
            integer, intent(in) :: istep                !! current step <br>
                                                        !! 当前步数
            real(rk), intent(in) :: t                   !! current time <br>
                                                        !! 当前时间
        end subroutine NS_fcn
    end interface

contains

    !> Navier-Stokes time integration <br>
    !> 在此实现所有的力学原理，必须通过力学方程合理计算加速度 (加速度由位移、速度、密度等决定，从而推导时间步迭代)
    subroutine navier_stokes_solver(region, istep, t)
        type(region_t), intent(inout) :: region     !! region of particle macroscopic properties <br>
                                                    !! 粒子宏观变量域
        integer, intent(in) :: istep                !! current time step <br>
                                                    !! 当前时间步
        real(rk), intent(in) :: t                   !! current time <br>
                                                    !! 当前时间
        type(list_t), allocatable :: grids(:, :)    !! 背景网格
        type(twin_list_t), save :: pairs            !! 粒子对，包含粒子对索引 item、核函数值 w、核函数导数值 dwdx

        ! ---------------------------- 初始化计算域内粒子的加速度 ------------------------- !
        region%particles%acc(:, :) = 0.0_rk                   ! 初始化加速度

        ! ------------------------------- 近邻粒子搜索（NNPS） --------------------------- !
        ! 近邻搜索，网格边长为 scale 倍光滑长度，可以使得中心网格中的粒子天然相互交互，相邻网格的粒子可能交互，且所有交互对在相邻网格内。
        ! 二维，单个网格内的粒子大致有 1 个 （+2 裕度）。@todo 改善预估粒子数量
        call build(region%particles%loc, scale_hsml=skf%scale_hsml, estimated_capacity=4, grid=grids)
        call query(region%particles%loc, grid=grids, scale=skf%scale, &
                   estimated_capacity=estimated_capacity, twin_list=pairs)

        !$omp parallel sections num_threads(cli%num_threads)
        !$omp section
        deallocate (grids)   ! 22.1% 耗时/ CPU 消耗

        !$omp section
        ! ------------------------------ 显示粒子配对、时间步信息 ------------------------- !
        if (mod(istep, nsave) == 0) write (99, "(*(a,g0.4))") "NOPP = ", pairs%log, &
            ", ANOP(+) = ", &
            real(pairs%log*2)/nums(4) + 1, ", t(s) = ", t, ", istep = ", istep
        ! @todo 进度条设计：在屏幕打印 100% 进度。

        ! ------------------------------------- 密度求和 -------------------------------- !
        call density(region%particles, pairs, nums(4))                          ! 对所有粒子进行密度求和

        ! ---------------------------------- 加速度（力）求解 ---------------------------- !
        call internal_force(region%particles, pairs, region%particles%acc, nums(4))       ! 求解内力：粘性力、压力
        call external_force(region%particles, pairs, region%particles%acc)                      ! 求解外力：重力、边界排斥力
        if (art_visc) call artificial_viscosity(region%particles, pairs, acc=region%particles%acc)   ! 求解人工粘性

        !$omp end parallel sections
        call update_virtual_particles(region%particles, region%virtual_particles, nums(2), nums(1))    ! 更新虚拟粒子

    end subroutine navier_stokes_solver

end module sph_navier_stokes_solver
